1.1. Chargement du tableau principal

On charge le fichier des pays en 2018

don <- read.table(file = "data/africa_pays_2018/data/africa_pays_2018.csv", # nom du fichier et chemin d'accès
                  sep = ";",                     # séparateur (ici, des points-virgule)
                  dec=",",                       # Type de décimale
                  header = TRUE,                 # ligne d'en-tête avec le nom des variables
                  encoding="UTF-8")              # encodage adapté au français

1.2 Choix des deux variables à analyser

En dehors de iso3 et nom, on ne garde que deux variables que l’on renomme X et Y avec colnames() et que l’on convertit en type numérique général. Il suffira par la suite de modifier le choix des variables X et Y pour faire d’autres analyses.

sel<-don[,c("iso3","nom","URBANI","INTERN")]
colnames(sel)<-c("CODE","NOM","X","Y")
head(sel)

  CODE                  NOM    X    Y
1  AGO               Angola 65.9 14.3
2  BDI              Burundi 13.2  2.7
3  BEN                Bénin 47.6 20.0
4  BFA         Burkina Faso 29.7 16.0
5  BWA             Botswana 69.8 47.0
6  CAF Rep. Centrafricaine  41.6  4.3

1.3 On est malin …

Mais comme on ne sait plus ce que sont X et Y, on le précise avec des chaînes de caractères qu’on pourra utiliser dans les graphiques. Et on peut préparer une version multilangue …

# Pour la version française
titre <- "Les pays d'Afrique en 2018"
nomX <- "Taux d'urbanisation (%)"
nomY <- "Taux d'accès à internet (%)"
auteur <- "Ecole d'été AfromapR, Bouaké 2025"

2.1 Vérification de la normalité de X et Y

La régression linéaire met en relation deux variables quantitatives X et Y dont on suppose que la distribution est normale (gaussienne) , c’est-à-dire unimodale et symérique.

2.1 Vérification de la normalité de X et Y

On peut tester la normalité des disributions par inspection visuelle à l’aide de hist()

On voit que X est assez symétrique et unimodale. En revanche Y est plutôt dissymétrique à gauche.

2.1 Vérification de la normalité de X et Y

Une analyse rapide avec boxplot() permet de repérer s’il y a des valeurs exceptionnelles.

par(mfrow=c(1,2),mar=c(4,2,2,2))
boxplot(sel$X, 
     col="lightyellow",
     horizontal = TRUE,
     main=nomX,
     cex.main=0.7,
     cex.axis=0.5)
boxplot(sel$Y, 
     col="lightyellow",
     horizontal = TRUE,
     main=nomY,
     cex.main=0.7, 
     cex.axis=0.5)

2.1 Vérification de la normalité de X et Y

Une analyse rapide avec boxplot() permet de repérer s’il y a des valeurs exceptionnelles.

Pour X comme Y, on ne voit pas apparaître de valeurs exceptionnelles c’est-à-dire de points situés au delà de l’intervalle défini par la boîte à moustache.

2.1 Vérification de la normalité de X et Y

Les fonctions qqnorm() et qqline() sont plus précises …

qqnorm(sel$X, col="blue",ylab=nomX)
qqline(sel$X,col="red")

X suit bien la droite et est proche d’une distribution gaussienne

2.1 Vérification de la normalité de X et Y

Les fonctions qqnorm() et qqline() sont plus précises …

qqnorm(sel$Y, col="blue",ylab=nomY)
qqline(sel$Y,col="red")

Y ne suit pas bien la droite donc elle n’est pas tout à fait gaussienne.

2.1 Vérification de la normalité de X et Y

Mais la solution la plus précise est le test de Shapiro qui pose l’hypothèse H0 : la distribution est normale.

shapiro.test(sel$X)

    Shapiro-Wilk normality test

data:  sel$X
W = 0.97546, p-value = 0.3931

shapiro.test(sel$Y)

    Shapiro-Wilk normality test

data:  sel$Y
W = 0.89916, p-value = 0.0005157

Si le test est supérieur à 0.05, on peut considérer que la distribution est gaussienne ce qui est le cas de X. Dans le cas contraire (p <0.05) la distribution n’est pas gaussienne ce qui est le cas de Y.

Mais en pratique on accepte en général de faire des analyses de régression avec des distributions non-gaussienne, surtout s’il n’y a pas de valeurs exceptionnelles.

2.2 Visualisation de la forme de la relation

On peut faire un simple plot(X,Y). Mais on peut aussi créer pour cela une fonction personalisée adapté à ses préférences

monplot <- function (varX , varY,  varN )
{ 
  Ymin <- min(varY)
  Ymax <- max(varY)
  Ysup <- Ymax+(Ymax-Ymin)*0.1
  plot(varX,varY,
     main = titre,      # titre
     cex.main = 1,      # police du titre
     cex = 0.6,         # taille des symboles
     pch = 19,          # cercles pleins
     ylim = c(Ymin, Ysup),
     col = "red")      # couleur des symboles
  text(varX,varY,varN,cex=0.5,pos=3) # nom des élément
  abline(v=mean(varX),lty=2,lwd=1,col="blue") # moyenne X
  abline(h=mean(varY),lty=2,lwd=1,col="blue") # moyenne Y
  grid()
  }

2.2 Visualisation de la forme de la relation

Je peux désormais utiliser ma fonction monplot() !

monplot(varX = sel$X,varY = sel$Y, varN = sel$CODE)

2.2 Visualisation de la forme de la relation

Je peux décider de ne pas afficher le label des points.

monplot(varX = sel$X,varY = sel$Y, varN = NULL)

2.3 Analyse de la corrélation

Je commence par calculer le coefficient de corrélation linéaire (r) et le pouvoir explicatif de X par rapport à Y (r2)

cor(sel$X,sel$Y)       # coefficient de corrélation (r)

[1] 0.530947

100*cor(sel$X,sel$Y)**2    # pouvoir explicatif (r2)

[1] 28.19048

2.3 Analyse de la corrélation

Puis, je teste la significativité de la corrélation linéaire …

cor.test(sel$X,sel$Y)  # test de significativité (p-value)

    Pearson's product-moment correlation

data:  sel$X and sel$Y
t = 4.2955, df = 47, p-value = 8.68e-05
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.2935825 0.7066417
sample estimates:
     cor 
0.530947 

2.3 Analyse de la corrélation

… et je la compare à celle du coefficient de corrélation de rang de Spearman

cor.test(sel$X,sel$Y, method="spearman")  # test de significativité (p-value)

    Spearman's rank correlation rho

data:  sel$X and sel$Y
S = 10474, p-value = 0.0007479
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.4656206 

Les deux coefficients sont assez proches (+0.53 pour Pearson et +0.47 pour Spearman) ce qui est rassurant. En général, si les deux coefficients diffèrent il y a un problème de valeur exceptionnelle ou de relation non linéaire.

2.3 Analyse de la corrélation

On peut conclure des analyses précédentes que :

• il existe une relation significative (p-value < 0.01)
• cette relation est positive (r = +0.53 )
• cette relation a un pouvoir explicatif assez faible (r2 = 28%)
• la relation est monotone et globalement linéaire car le coefficient de Spearman est à peu près égale au coefficient de Pearson. Il n’y a donc pas besoin de transformer les variables ou retirer des valeurs exceptionnelles.

3.1 Hypothèses statistiques

Conditions a priori

1. X et Y sont deux variables normales (gaussienne)
2. il existe une corrélation significative entre X et Y (p< 0.05)
3. X explique une part suffisamment forte de Y (r2 > 20% )
4. Le nuage de point affiche une forme linéaire
5. les points sont répartis de façon régulière le long du nuage de points
6. Il n’y a pas de valeurs exceptionnelles susceptibles de perturber le calcul.

3.1 Hypothèses statistiques

Méthode des moindres carrés ordinaire (MCO)

• La droite \(y_i = a.x_i + b + \epsilon_i\) qui minimise la somme des carrés des écarts entre les valeurs observées \(y_i\) et les valeurs estimées \(\hat{y_i}\) a pour équation :

• \(COV(X,Y) = \sum_{i=1}^k \sum_{j=1}^k (x_{i}-\bar{x})^2.(y_{i}-\bar{y})^2\)

• \(a = COV(X,Y) / (\sigma_X)^2\)

• \(b = \bar{y} - a.\bar{x}\)

3.1 Hypothèses statistiques

Analyse de la variance

• La somme des carré des écarts totale (\(SCE_{tot}\)) correspond à la variance de la variable à expliquer : \(SCE_{tot} = \sum_{i=1}^k (y_{i}-\bar{y})^2\)

• La somme des carré des écarts résiduels (\(SCE_{err}\)) correspond à la somme des carrés des différences entre valeurs observées et estimées \(SCE_{error} = \sum_{i=1}^k (y_{i}-\hat{y})^2\)

• Le pouvoir explicatif d’un modèle de régression correspond à la part de la variance de Y expliquée par X.

• \(Var. expliquée = (SCE_{tot}-SCE_{res}) / SCE_{tot} = r(X,Y)^{2}\)

3.1 Hypothèses statistiques

Vérifications a posteriori

Un modèle de régression n’est valide que si les résidus de ce modèle \(\epsilon_i = (y_i - \hat{y}_i)\) remplissent les conditions suivantes :

1. Normalité de la distribution des résidus
2. Absence d’autocorrélation des résidus
3. Homogénéité de la variance des résidus
4. Absence de valeur à fort effet de levier

Si ces quatre conditions ne sont pas remplies, les estimations de Y en fonction de X seront entâchées d’erreur et leur intervalle de confiance ne sera pas valable.

3.2 La fonction lm()

La fonction lm() ou lm est l’abbréviation de linear model permet d’effectuer la plupart des modèles de régression linéaire basés sur la méthode des moindres carrés ordinaire. Sa syntaxe est a priori très simple et renvoie les coefficients b et a du modèle de régression.

lm(sel$Y~sel$X)

Call:
lm(formula = sel$Y ~ sel$X)

Coefficients:
(Intercept)        sel$X  
     1.1223       0.5344  

3.2 La fonction lm()

Mais en réalité lm() crée une liste de résultats que l’on a intérêt à stocker pour en examiner les composantes une à une.

monmodel<-lm(sel$Y~sel$X)
str(monmodel)

List of 12
 $ coefficients : Named num [1:2] 1.122 0.534
  ..- attr(*, "names")= chr [1:2] "(Intercept)" "sel$X"
 $ residuals    : Named num [1:49] -22.042 -5.477 -6.562 -0.995 8.574 ...
  ..- attr(*, "names")= chr [1:49] "1" "2" "3" "4" ...
 $ effects      : Named num [1:49] -177.69 69.57 -3.8 2.06 10.97 ...
  ..- attr(*, "names")= chr [1:49] "(Intercept)" "sel$X" "" "" ...
 $ rank         : int 2
 $ fitted.values: Named num [1:49] 36.34 8.18 26.56 17 38.43 ...
  ..- attr(*, "names")= chr [1:49] "1" "2" "3" "4" ...
 $ assign       : int [1:2] 0 1
 $ qr           :List of 5
  ..$ qr   : num [1:49, 1:2] -7 0.143 0.143 0.143 0.143 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:49] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:2] "(Intercept)" "sel$X"
  .. ..- attr(*, "assign")= int [1:2] 0 1
  ..$ qraux: num [1:2] 1.14 1.27
  ..$ pivot: int [1:2] 1 2
  ..$ tol  : num 1e-07
  ..$ rank : int 2
  ..- attr(*, "class")= chr "qr"
 $ df.residual  : int 47
 $ xlevels      : Named list()
 $ call         : language lm(formula = sel$Y ~ sel$X)
 $ terms        :Classes 'terms', 'formula'  language sel$Y ~ sel$X
  .. ..- attr(*, "variables")= language list(sel$Y, sel$X)
  .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:2] "sel$Y" "sel$X"
  .. .. .. ..$ : chr "sel$X"
  .. ..- attr(*, "term.labels")= chr "sel$X"
  .. ..- attr(*, "order")= int 1
  .. ..- attr(*, "intercept")= int 1
  .. ..- attr(*, "response")= int 1
  .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..- attr(*, "predvars")= language list(sel$Y, sel$X)
  .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. ..- attr(*, "names")= chr [1:2] "sel$Y" "sel$X"
 $ model        :'data.frame':  49 obs. of  2 variables:
  ..$ sel$Y: num [1:49] 14.3 2.7 20 16 47 4.3 46.8 23.2 8.6 8.7 ...
  ..$ sel$X: num [1:49] 65.9 13.2 47.6 29.7 69.8 41.6 51 56.7 44.8 67.2 ...
  ..- attr(*, "terms")=Classes 'terms', 'formula'  language sel$Y ~ sel$X
  .. .. ..- attr(*, "variables")= language list(sel$Y, sel$X)
  .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:2] "sel$Y" "sel$X"
  .. .. .. .. ..$ : chr "sel$X"
  .. .. ..- attr(*, "term.labels")= chr "sel$X"
  .. .. ..- attr(*, "order")= int 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 1
  .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. ..- attr(*, "predvars")= language list(sel$Y, sel$X)
  .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. .. ..- attr(*, "names")= chr [1:2] "sel$Y" "sel$X"
 - attr(*, "class")= chr "lm"

3.3 Analyse du modèle

Un résumé des résultats principaux est fourni avec summary() appliqué à l’objet créé par lm().

summary(monmodel)

Call:
lm(formula = sel$Y ~ sel$X)

Residuals:
    Min      1Q  Median      3Q     Max 
-28.337 -11.417  -2.522  12.700  33.105 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.1223     6.1037   0.184    0.855    
sel$X         0.5344     0.1244   4.295 8.68e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 16.2 on 47 degrees of freedom
Multiple R-squared:  0.2819,    Adjusted R-squared:  0.2666 
F-statistic: 18.45 on 1 and 47 DF,  p-value: 8.68e-05

On obtient ainsi :

• l’équation de la droite Y = a.X+b
• la significativité et l’intervalle de confiance de a et b
• le pouvoir explicatif du modèle \(r(X,Y)^2\)

3.4 Tracé de la droite de régression

On peut tracer la droite de régression avec abline()

monplot(sel$X,sel$Y,sel$CODE)
abline(monmodel, col="blue",lwd=2)

Cette droite a pour équation Y = 0.5X + 2.8 ce qui veut dire que chaque fois que l’urbanisation augemente de 1 pt l’usage d’internet augmente de 0.5 pts. Un pays qui aurait une urbanisation de 50% devrait donc avoir un taux d’usage d’internet égal à 0.5 x 50 + 2.8 = 27.8%

3.5 Analyse des résidus

On peut extraire de l’objet créé par lm() les valeurs estimées de Y et les résidus c’est-à-dire les erreurs d’estimation.

sel$Y_estim<-monmodel$fitted.values # Valeurs estimées
sel$Y_resid<-monmodel$residuals     # Valeurs résiduelles
tri <- sel[order(sel$Y_resid),]     # Tri du tableau en fonction des résidus
tri

   CODE                  NOM    X    Y   Y_estim     Y_resid
10  COG                Congo 67.2  8.7 37.036626 -28.3366260
40  SOM              Somalie 45.3  2.0 25.332411 -23.3324109
24  LBY                Libye 80.3 21.8 44.037777 -22.2377775
1   AGO               Angola 65.9 14.3 36.341855 -22.0418552
14  ERI             Erythrée 40.4  1.3 22.713660 -21.4136596
23  LBR              Libéria 51.4  8.0 28.592489 -20.5924891
20  GNB        Guinée-Bissau 43.6  3.9 24.423865 -20.5238646
6   CAF Rep. Centrafricaine  41.6  4.3 23.354986 -19.0549865
9   COD  Congo, Rép. dém. du 44.8  8.6 25.065191 -16.4651914
39  SLE         Sierra Leone 42.3  9.0 23.729094 -14.7290938
19  GMB               Gambie 61.6 19.8 34.043767 -14.2437674
21  GNQ   Guinée équatoriale 72.4 26.2 39.815709 -13.6157090
27  MDG           Madagascar 37.6  9.8 21.217230 -11.4172303
43  TGO                 Togo 42.0 12.4 23.568762 -11.1687621
28  MLI                 Mali 42.8 13.0 23.996313 -10.9963133
29  MOZ           Mozambique 36.3 10.0 20.522460 -10.5224596
48  ZMB               Zambie 43.8 14.3 24.530752 -10.2307524
30  MRT           Mauritanie 54.1 20.8 30.035475  -9.2354745
8   CMR             Cameroun 56.7 23.2 31.425016  -8.2250160
42  TCD                Tchad 23.2  6.5 13.521308  -7.0213081
3   BEN                Bénin 47.6 20.0 26.561621  -6.5616207
2   BDI              Burundi 13.2  2.7  8.176918  -5.4769176
33  NER                Niger 16.5  5.3  9.940566  -4.6405665
37  SSD           Sud Soudan 19.8  8.0 11.704215  -3.7042153
18  GIN               Guinée 36.3 18.0 20.522460  -2.5224596
4   BFA         Burkina Faso 29.7 16.0 16.995162  -0.9951619
22  KEN                Kenya 27.3 17.8 15.712508   2.0874918
31  MWI               Malawi 17.1 13.8 10.261230   3.5387701
45  TZA             Tanzanie 34.2 25.0 19.400138   5.5998624
15  ETH             Ethiopie 21.0 18.6 12.345542   6.2544578
17  GHA                Ghana 56.4 39.0 31.264684   7.7353157
5   BWA             Botswana 69.8 47.0 38.426168   8.5738325
49  ZWE             Zimbabwe 32.2 27.1 18.331259   8.7687405
46  UGA              Ouganda 24.1 23.7 14.002303   9.6976968
36  SDN               Soudan 34.8 30.9 19.720801  11.1791990
35  RWA               Rwanda 17.3 21.8 10.368118  11.4318823
25  LSO              Lesotho 28.4 29.0 16.300391  12.6996089
11  DJI             Djibouti 77.9 55.7 42.755124  12.9448762
16  GAB                Gabon 89.6 62.0 49.008061  12.9919394
34  NGA              Nigéria 50.8 42.0 28.271826  13.7281743
7   CIV        Côte d'Ivoire 51.0 46.8 28.378713  18.4212865
47  ZAF       Afrique du Sud 66.7 56.2 36.769406  19.4305935
38  SEN              Sénégal 47.5 46.0 26.508177  19.4918232
12  DZA              Algérie 72.9 59.6 40.082929  19.5170714
32  NAM              Namibie 50.5 51.0 28.111494  22.8885060
13  EGY               Egypte 42.7 46.9 23.942869  22.9571306
44  TUN              Tunisie 69.1 64.2 38.052060  26.1479398
26  MAR                Maroc 62.8 64.8 34.685094  30.1149058
41  SWZ            Swaziland 23.9 47.0 13.895415  33.1045846

3.5 Analyse des résidus

Les résidus négatifs correspondent aux pays qui ont moins accès à internet que ce que leur urbanisation laisserait prévoir :

head(tri,5)

   CODE      NOM    X    Y  Y_estim   Y_resid
10  COG    Congo 67.2  8.7 37.03663 -28.33663
40  SOM  Somalie 45.3  2.0 25.33241 -23.33241
24  LBY    Libye 80.3 21.8 44.03778 -22.23778
1   AGO   Angola 65.9 14.3 36.34186 -22.04186
14  ERI Erythrée 40.4  1.3 22.71366 -21.41366

• Ainsi le Congo n’a que 8.7% d’accès à internet (Y) alors que son taux d’urbanisation est de 67.2% (X) ce qui laisserait prévoir un taux d’accès à internet de 37.0% (Y_estim). Il a donc un résidus de -28 pts (Y_resid) qui signifie que son taux d’accès à internet est beaucoup plus faible que ce que laisserait prévoir son urbanisation.

3.5 Analyse des résidus

Les résidus positifs correspondent aux pays qui ont plus accès à internet que ce que leur urbanisation laisserait prévoir :

tail(tri,5)

   CODE       NOM    X    Y  Y_estim  Y_resid
32  NAM   Namibie 50.5 51.0 28.11149 22.88851
13  EGY    Egypte 42.7 46.9 23.94287 22.95713
44  TUN   Tunisie 69.1 64.2 38.05206 26.14794
26  MAR     Maroc 62.8 64.8 34.68509 30.11491
41  SWZ Swaziland 23.9 47.0 13.89542 33.10458

-Ainsi le Swaziland a 47 % d’accès à internet (Y) alors que son taux d’urbanisation est de 23.9% (X) ce qui laisserait prévoir un taux d’accès à internet de 13.9% (Y_estim). Il a donc un résidus de +33 pts (Y_resid) qui signifie que son taux d’accès à internet est beaucoup plus fort que ce que laisserait prévoir son urbanisation.

4.1 Hypothèses statistiques

Vérifications a posteriori

Un modèle de régression n’est valide que si les résidus de ce modèle \(\epsilon_i = (y_i - \hat{y}_i)\) remplissent les conditions suivantes :

1. Normalité de la distribution des résidus
2. Absence d’autocorrélation des résidus
3. Homogénéité de la variance des résidus
4. Absence de valeur à fort effet de levier

Si ces quatre conditions ne sont pas remplies, les estimations de Y en fonction de X seront entâchées d’erreur et leur intervalle de confiance ne sera pas valable.

On charge le package car (companion to applied regession) pour réaliser certains des diagnostics

library(car)

4.1 Diagnostic 1 : Indépendance des résidus ?

L’objectif est de savoir si les résidus se répartissent régulièrement de part et d’autre de la droite de régression tout au long de celle-ci. Si c’est bien le cas le graphique residuals Vs Fitted généré par plot(monmodel,1) devrait donner une droite horizontale :

plot(monmodel,1,labels.id = sel$CODE)

4.1 Diagnostic 1 : Indépendance des résidus ?

4.1 Diagnostic 1 : Indépendance des résidus ?

On peut tester statistiquement l’indépendance des résidus à l’aide du test de Durbin-Watson qui mesure si deux valeurs successives ont des résidus proches. L’indépendance des résidus est rejetée si p-value < 0.05

durbinWatsonTest(monmodel)

 lag Autocorrelation D-W Statistic p-value
   1      0.00312332       1.94811    0.83
 Alternative hypothesis: rho != 0

Ici on trouve p-value > 0.05 donc les résidus sont indépendants.

Diagnostic 2 : Normalité des résidus ?

L’objectif est de savoir si les résidus ont une distribution normale Si c’est bien le cas le graphique généré par plot(monmodel,2) devrait donner une droite oblique :

plot(monmodel,2,labels.id = sel$CODE)

Diagnostic 2 : Normalité des résidus ?

Diagnostic 2 : Normalité des résidus ?

On peut tester statistiquement la normalité des résidus à l’aide du test de Shapiro-Wilk. Les résidus sont normaux si p-value > 0.05

shapiro.test(monmodel$residuals)

    Shapiro-Wilk normality test

data:  monmodel$residuals
W = 0.96515, p-value = 0.1543

Ici on trouve une p-value supérieure à 0.05 donc la distribution des résidus est approximativement gaussienne.

Diagnostic 3 : Homogénéité des résidus ?

L’objectif est de savoir si la variance des résidus est constante, c’est-à-dire si il s’écarte environ de la même distance tout au long de la droite . Si c’est bien le cas le graphique généré par plot(monmodel,3) devrait donner une droite horizontale

plot(monmodel,3,labels.id = sel$CODE)

Diagnostic 3 : Homogénéité des résidus ?

Diagnostic 3 : Homogénéité des résidus ?

On peut tester statistiquement l’homogénéité des résidus à l’aide du test de Breush-Pagan. L’hypothèse d’homogénéité est rejetée si la p-value est inférieure à 0.05.

ncvTest(monmodel)

Non-constant Variance Score Test 
Variance formula: ~ fitted.values 
Chisquare = 3.029338, Df = 1, p = 0.081771

Ici, la p-value est supérieure à 0.05 donc les résidus sont globalement homogènes.

Diagnostic 4 : Absence de valeur exceptionnelles ?

L’objectif est de savoir s’il existe des valeurs qui exercent une influence exceptionnelle sur les résultats de la régression. On peut reprérer ces valeurs de plusieurs manières, notamment à l’aide de la distance de Cook générée par plot(monmodel,4).

plot(monmodel,4,labels.id = sel$CODE)

Diagnostic 4 : Absence de valeur exceptionnelles ?

Les pays qui influent sur la régression sont ceux qui sont à la fois éloignés de la droite et situés aux extrémités de celle-ci. Si on les retirait, la pente de la droite changerait.C’est ici le cas du Congo, de la Libye et du Swaziland

Diagnostic 4 : Absence de valeur exceptionnelles ?

Le test statistique de Bonferroni permet de déterminer si les pays à fort effet de levier sont réellement de nature à pertureber les résultats.

outlierTest(monmodel, labels = sel$CODE)

No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:
    rstudent unadjusted p-value Bonferroni p
SWZ 2.176162           0.034715           NA

Ici, on voit que le Swaziland est le pays le plus influent. Mai le test de Bonferoni renvoie NA et il n’est donc pas considéré comme réellement exceptionnel.

Bilan des diagnostics :

En résumé, nous avons pu vérifier que notre modèle était correct sur chacun des quatre diagnostics :

1. Absence d’autocorrélation des résidus
2. Normalité des résidus
3. Homogénéité des résidus
4. Absence de valeur exceptionnellement influente.

Visualisation rapide

Lorsqu’on utilise le package car, on peut visualiser rapidement une relation en utilisant un petit nombre d’instruction avec la fonction scatterplot() du package car. Le paramètre smooth=TRUE permet de visualiser la dispersion des points sous forme d’ellipse :

scatterplot(formula = sel$X~sel$Y, 
            smooth=FALSE,
            ellipse=TRUE,
            xlab=nomX,
            ylab=nomY,
            main = "Affichage de l'ellipse de dispersion")

Visualisation rapide

Lorsqu’on utilise le package car, on peut visualiser rapidement une relation en utilisant un petit nombre d’instruction avec la fonction scatterplot() du package car. Le paramètre ellipse=TRUE permet de visualiser la dispersion des points sous forme d’un ajustement non linéaire par la méthode loess.

scatterplot(formula = sel$X~sel$Y, 
            smooth=T,
            ellipse=F,
            xlab=nomX,
            ylab=nomY,
            main = "Affichage de l'ajustement non linéaire (loess)")